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Abstract 

The idea of stable, localized bundles of energy has strong appeal as a model for particles. In 
the 1950s, John Wheeler envisioned such bundles as smooth configurations of electromagnetic 
energy that he called geons, but none were found. Instead, particle-like solutions were found 
in the late 1960s with the addition of a scalar field, and these were given the name boson stars. 
Since then, boson stars find use in a wide variety of models as sources of dark matter, as black 
hole mimickers, in simple models of binary systems, and as a tool in finding black holes in 
higher dimensions with only a single killing vector. We discuss important varieties of boson 
stars, their dynamic properties, and some of their uses, concentrating on recent efforts. 
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1 Introduction 



Particle-like objects have a very long and broad history in science, arising long before Newton's 
corpuscles of light, and spanning the range from fundamental to astronomical. In the mid-1950s, 
John Wheeler sought to construct stable, particle-like solutions from only the smooth, classical 
fields of electromagnetism coupled to general relativity [2211 1183] . Such solutions would represent 
something of a "gravitational atom," but the solutions Wheeler found, which he called geons, were 
unstable. However in the following decade, Kaup replaced electromagnetism with a complex scalar 
field [128] . and found Klein-Gordon geons that, in all their guises, have become well-known as 
today's boson stars (see Section II of [195] for a discussion of the naming history of boson stars). 

As compact, stationary configurations of scalar field bound by gravity, boson stars are called 
upon to fill a number of different roles. Most obviously, could such solutions actually represent 
astrophysical objects, either observed directly or indirectly through its gravity? Instead, if con- 
structed larger than a galaxy, could a boson star serve as the dark matter halo that explains the 
flat rotation curve observed for most galaxies? 

The equations describing boson stars are relatively simple, and so even if they do not exist in 
nature, they still serve as a simple and important model for compact objects, ranging from particles 
to stars and galaxies. In all these cases, boson stars represent a balance between the dispersive 
nature of the scalar field and the attraction of gravity holding it together. 

This review is organized as follows. The rest of this section describes some general features 
about boson stars. The system of equations describing the evolution of the scalar field and grav- 
ity (i.e., the Einstein-Klein-Gordon equations) are presented in Section [2] These equations are 
restricted to the spherical symmetric case (with a harmonic ansatz for the complex scalar field 
and a simple massive potential) to obtain a boson star family of solutions. To accommodate all 
their possible uses, a large variety of boson star types have come into existence, many of which 
are described in more detail in Section [3] For example, one can vary the form of the scalar field 
potential to achieve a large range of masses and compactness than with just a mass term in the 
potential. Certain types of potential admit soliton-like solutions even in the absence of gravity, 
leading to so-called Q-stars. One can adopt Newtonian gravity instead of general relativity, or con- 
struct solutions from a real scalar field instead of a complex one. It is also possible to find solutions 
coupled to an electromagnetic field or a perfect fluid, leading respectively to charged boson stars 
and fermion-boson stars. Rotating boson stars are found to have an angular momentum which is 
not arbitrary, but instead quantized. Multi-state boson stars with more than one complex scalar 
field are also considered. 

We discuss the dynamics of boson stars in Section |U Arguably the most important property of 
boson star dynamics concerns their stability. A stability analysis of the solutions can be performed 
either by studying linear perturbations, catastrophe theory, or numerical non-linear evolutions. 
The latter option allows for the study of the final state of perturbed stars. Possible endstates 
include dispersion to infinity of the scalar field, migration from unstable to stable configurations, 
and collapse to a black hole. There is also the question of formation of boson stars. Full numerical 
evolutions in 3D allow for the merger of binary boson stars, which display a large range of different 
behaviours as well producing distinct gravitational wave signatures. 

Finally, we review the impact of boson stars in astronomy in Section [5] (as an astrophysical 
object, black hole mimickers and origin of dark matter) and in mathematics in Section [6] (studies of 
critical behaviour, the Hoop conjecture and higher dimensions). We conclude with some remarks 
and future directions. 
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1.1 The nature of a boson star 

Boson stars are constructed with a complex scalar field coupled to gravity (as described in Sec- 
tion [2]). A complex scalar field </>(t,r) can be decomposed into two real scalar fields </>r and 4>\ 
mapping every spacetime event to the complex plane 

4>{t,T) = <f> K {t,T)+i<t>i{t,r). (1) 

Such a field possesses energy because of its spatial gradients and time derivatives and this energy 
gravitates holding the star together. Less clear is what supports the star against the force of gravity. 
Its constituent scalar field obeys a Klein-Gordon wave equation which tends to disperse fields. 
This is the same dispersion which underlies the Heisenberg Uncertainty Principle. Indeed, Kaup's 
original work |128j found energy eigenstates for a semi-classical, complex scalar field, discovering 
that gravitational collapse was not inevitable. Ruffini and Bonazzola [189] followed up on this work 
by quantizing a real scalar field representing some number of bosons and they found the same field 
equations. 

None of this guarantees that such solutions balancing dispersion against gravitational attraction 
exist. In fact, a widely known theorem, Derrick's Theorem [69] (see also |187j ). uses a clever 
scaling argument to show that no regular, static, nontopological localized scalar field solutions are 
stable in three (spatial) dimensional flat space. This constraint is avoided by adopting a harmonic 
ansatz for the complex scalar field 

0(r,t) = ^ o (r)e <ut (2) 

and by working with gravity. Although the field is no longer static, as shown in Section [2] the 
spacetime remains static. The star itself is a stationary, soliton-like solution as demonstrated in 
Figure [1] 

There are, of course, many other soliton and soliton-like solutions in three dimensions finding 
a variety of ways to evade Derrick's Theorem. For example, the field-theory monopole of 't Hooft 
and Polyakov is a localized solution of a properly gauged triplet scalar field. Such a solution is 
a topological soliton because the monopole possesses false vacuum energy which is topologically 
trapped. The monopole is one among a number of different topological defects that requires an 
infinite amount of energy to "unwind" the potential energy trapped within (see [219] for a general 
introduction to defects and the introduction of |190j for a discussion of relevant classical field theory 
concepts) . 

In Section [21 we present the underlying equations and mathematical solutions, but here we are 
concerned with the physical nature of these boson stars. When searching for an actual boson star, 
we look not for a quantized wave function, or even a semiclassical one. Instead, we search for a 
fundamental scalar, say the long-sought Higgs boson. The Large Hadron Collider (LHC) hopes to 
determine the existence and nature of the Higgs, with evidence at the time of writing suggesting 
a Higgs boson with mass ~ 125 GeV/c 2 [185] . If the Higgs does not ultimately appear, there 
are other candidates such as an axion particle. Boson stars are then either a collection of stable 
fundamental bosonic particles bound by gravity, or else a collection of unstable particles that, with 
the gravitational binding, have an inverse process efficient enough to reach an equilibrium. They 
can thus be considered a Bose-Einstein condensate (BEC), although boson stars can also exist in 
an excited state as well. 

Indeed, applying the uncertainty principle to a boson star by assuming it to be a macroscopic 
quantum state results in an excellent estimate for the maximum mass of a BS. One begins with 
the Heisenberg uncertainty principle of quantum mechanics 

Ap Ax > h (3) 
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and assumes the BS is confined within some radius Ax = R with a maximum momentum of 
Ap = mc where m is the mass of the constituent particle 



mcR > h. (4) 

This inequality is consistent with the star being described by a Compton wavelength of Ac = 
h/ {mc). We look for the maximum possible mass M max for the boson star which will saturate the 
uncertainty bound and drive the radius of the star towards its Schwarzschild radius Rg = 2GM/c 2 . 
Substituting yields 

2GmM max 

= K (5) 

c 

which gives an expression for the maximum mass 

1 he 

M ^=2G^- (6) 



Recognizing the Planck mass Mpi an ck = ^hc/G, we obtain the estimate of M max = 0.5Mp lanck /m. 
This simple estimate indicates that the maximum mass of the BS is inversely related to the mass of 
the constituent scalar field. We will see below in Section[2]that this inverse relationship continues to 
hold with the explicit solution of the differential equations for a simple mass term in the potential, 
but can vary with the addition of self-interaction terms. Indeed depending on the strength of the 
coupling m and the other parameters of the self-interaction potential, the size and mass of the 
boson stars can vary from atomic to astrophysical scales. 

Despite their connection to fundamental physics, one can also view boson stars in analogy with 
models of neutron stars. In particular, as we discuss in the following sections, both types of stars 
demonstrate somewhat similar mass versus radius curves for their solutions with a transition in 
stability at local maxima of the mass. There is also a correspondence between (massless) scalar 
fields and a stiff, perfect fluid (see Section 2.1 and Appendix A of Ref. [36]), but the correspondence 
does not mean that the two are equivalent [79]. More than just an analogy, boson stars can serve as 
a very useful model of a compact star, having certain advantages over a fluid neutron star model: 
(i) the equations governing its dynamics avoid developing discontinuities, in particular there is 
no sharp stellar surface, (ii) there is no concern about resolving turbulence, and (iii) one avoids 
uncertainties in the equation of state. 



1.2 Other reviews 

A number of other reviews of boson stars have appeared. Most recently, Schunck and Miclkc [195 
concentrate on the possibility of detecting BS, extending their previous reviews [163, 164 . In 1992, 
a number of reviews appeared: Jetzer |124] concentrates on the astrophysical relevance of BS (in 
particular their relevance for explaining dark matter) while Liddle and Madsen [154] focus on their 
formation. Other reviews include |209[ 1151] . 
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Figure 1: Demonstration of the solitonic nature of the (mini-)boson star. Shown are snapshots 
of the magnitude squared of the complex scalar held for a head-on collision of two identical mini- 
boson stars. The interacting stars display an interference pattern as they pass through each other, 
recovering their individual identities after the collision. However, note that the BSs have a larger 
amplitude after their interaction and so are not true solitons. The collision can therefore be 
considered inelastic. From [50]. See also [143] (e.g., Figure 5.12). 
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2 Solving for Boson Stars 



In this section, we present the equations governing boson star solutions, namely the Einstein 
equations for the geometry description and the Klein-Gordon equation to represent the (complex) 
scalar field. We refer to this coupled system as the Einstein-Klein-Gordon (EKG) equations. 

The covariant equations describing boson stars are presented in Section [2T21 which is followed 
by choosing particular coordinates consistent with a 3+1 decomposition in Section [2. 31 A form for 
the potential of the scalar field is then chosen and solutions are presented in Section 12.41 



2.1 Conventions 

Throughout this review, Roman letters from the beginning of the alphabet a, b, c, . . . denote space- 
time indices ranging from to 3, while letters near the middle i, j,k, . . . range from 1 to 3, denoting 
spatial indices. Unless otherwise stated, we use units such that h = c = 1 so that the Planck mass 
becomes Mpi anc k = G -1 / 2 . We also use the signature convention (— , +, +, +) for the metric. 

2.2 The Lagrangian, evolution equations and conserved quantities 

The EKG evolution equations can be derived from the action [220 

S = I (t^ R + £m) V~9^x (7) 



16nG 

where R is the Ricci scalar of the spacetime represented by the metric g a b, and its determinant 
\J—g. The term Cm describes the matter which here is that of a complex scalar field, 4> 



r 1 r 



g ah V a 4>V b (j> + V (|0| 2 )] (8) 



where 4> is the complex conjugate of the field and y(|<^)| 2 ) a potential depending only on the 
magnitude of the scalar field, consistent with the U(l) invariance of the field in the complex plane. 

Variation of the action in Eq. (J7]) with respect to the metric g ab leads to the well-known Einstein 
equations 

R 

Rab - -^9ab = 8TrGT ab (9) 

Tab = \ [V a 4>Vb4> + V a <}>Vb4>] -\9ab[9 Cd Vc4>Vd<t>+V(\ct>\ 2 )] , (10) 

where R ao is the Ricci tensor and T a b is the real stress-energy tensor. Eqs. (J9j> form a system of 
10 non-linear partial differential equations for the spacetime metric components g a b coupled to the 
scalar field via the stress-energy tensor given in Eq. (|10j) . 

On the other hand, the variation of the action in Eq. ([7|) with respect to the scalar field <j>, leads 
to the Klein-Gordon (KG) equation 

An equivalent equation is obtained when varying the action with respect to the complex conjugate 
(p. The simplest potential leading to boson stars is the so-called free field case, where the potential 
takes the form 

V{m=m 2 H 2 , (12) 
with m a parameter that can be identified with the bare mass of the field theory. 
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According to Noether's theorem, the invariance of the Klein-Gordon Lagrangian in Eq. (JSj> 
under global U(l) transformations <j> — > 4>e lv implies the existence of a conserved current 

Ja = ^(^a<p-<f>^a4>), (13) 

satisfying the conservation law 

V a J a = -^d a (^g~g ab J b ) =0. (14) 

The spatial integral of the time component of this current defines the conserved Noether charge, 
given by 

N = J g 0b J b V=gdx 3 , (15) 
which can be associated with the total number of bosonic particles [189 . 

2.3 The 3+1 decomposition of the spacetime 

Although the spacetime description of general relativity is very elegant, the covariant form of 
Einstein equations is not suitable to describe how an initial configuration evolves towards the 
future. It is therefore more intuitive to instead consider a succession of spacetime geometries, 
where the evolution of a given slice is given by the Einstein equations (for more detailed treatments 
see [5] [25j [35] 98 ) . In order to convert the four-dimensional, covariant Einstein equations to a 
more intuitive "space+time" or 3+1 decomposition, the following steps are taken: 

• specify the choice of coordinates. The spacetime is foliated by a family of spacelike hyper- 
surfaces, which are crossed by a congruence of time lines that will determine our observers 
(i.e., coordinates). This congruence is described by the vector field t a = an a + (3 a , where a 
is the lapse function which measures the proper time of the observers, f3 a is the shift vector 
that measures the displacement of the observers between consecutive hypersurfaces and n a 
is the timelike unit vector normal to the spacelike hypersurfaces. 

• decompose every J^D object into its 3+1 components. The choice of coordinates allows for 
the definition of a projection tensor j a b = 8% + n a n b . Any four-dimensional tensor can be 
decomposed into 3+1 pieces using the spatial projector to obtain the spatial components, or 
contracting with n a for the time components. For instance, the line element can be written 
in a general form as 

ds 2 = -a 2 dt 2 + j ij (dx i + ftdt) (dx j + fiHt) . (16) 
The stress-energy tensor can then be decomposed into its various components as 

T = T ab n a n b , Si=T ab n a Yi , =T abl a n b y (17) 

• write down the field equations in terms of the 3+1 components. Within the framework 
outlined here, the induced (or equivalently, the spatial 3D) metric 7^ and the scalar field 
<f> are as yet still unknown (remember that the lapse and the shift just describe our choice 
of coordinates). In the original 3+1 decomposition (ADM formulation [TU]) an additional 
geometrical tensor Kij = — (l/2)£ n 7y = — 1/ (2a) (dt — is introduced to describe 
the change of the induced metric along the congruence of observers. Loosely speaking, one 
can view the determination of 7^ and Kij as akin to the specification of a position and 
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velocity for projectile motion. In terms of the extrinsic curvature and its trace, trK = K, 
the Einstein equations can be written as 



(trK) 2 - Ki j Kf = 16 7T G t 



Vj (Ki 3 - trK 8^) =8nGSi 



(18) 
(19) 

) 

(20) 

In a similar fashion, one can introduce a quantity Q = —£ n 4> for the Klein-Gordon equation 
which reduces it to an equation first order in time, second order in space 



(8 t - Cp) = -V<V 3 -a + a - 2K l k K jk + trKK K 



8ttG 



Si 



Hi 



(trS 



a V7 



(21) 



• enforce any assumed symmetries. Although the boson star is found by a harmonic ansatz 
for the time dependence, here we choose to retain the full time-dependence. However, a con- 
siderably simplification is provided by assuming that the spacetime is spherically symmetric. 
Following [143 , the most general metric in this case can be written in terms of spherical 
coordinates as 



ds 2 = (-a 2 + a 2 /3 2 ) dt 2 + 2 a 2 (3 dt dr + a 2 dr 2 + r 2 b 2 d£l 



(22) 



where a{t, r) is the lapse function, fi(t,r) is the radial component of the shift vector and 
a(t, r), b(t, r) represent components of the spatial metric, with dft 2 the metric of a unit two- 
sphere. With this metric, the extrinsic curvature only has two independent components 
Kj = diag (K r r ,K e g,K e g). The constraint equations, Eqs. (fTS]) and ([19]). can now be 
written as 



2 
arb 



d r (rb) 



1 

rb 



d r K e 



rb 

dr [ — d r (rb) 
1 a 



AK r r K 
8nG 



d r (rb) 
rb 



[K e e -K%) 



a z 
2nG 



2K\K% (23) 

|$| 2 + |n| 2 + a V(|0| 2 )(]24) 
+ ni>) (25) 



where we have defined the auxiliary scalar-field variables 

a 



$ =8 r 



n : 



(cW - [3d r 4>) 



The evolution equations for the metric and extrinsic curvature components reduce to 
dta = d r (a/3) — aaK r r 



8 t b = -d r (rb) - abK 6 
r 

d t K\ - (3d r K r r 
d t K\ - f3d r K e g 



-Or 



a 
a 



d r a 
a 



1 



(rb) 2 a(rb) 



arb 



-0 r 



d r (rb) 



d r (rb) 



trK K r 



4nG 



[2|$| 



a(trKK%-4:irGV(\(/>\ 2 )] 



(26) 

(27) 
(28) 

+ «MH 2 )]} 
(29) 
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Similarly, the reduction of the Klein-Gordon equation to first order in time and space leads 
to the following set of evolution equations 



-n 



d t $ = d r [p$ + -n 



1 



(rb) 



{rbf (pn + - 

V a 



P 



d r (rb) 



rb 



II — aa 



dV 

W 2 



(30) 
(31) 

(32) 



This set of equations, Eqs. (|23|) - (|32 |) . describes general, time-dependent, spherically sym- 
metric solutions of a gravitationally-coupled complex scalar field. In the next section, we 
proceed to solve for the specific case of a boson star. 



2.4 Mini-boson stars 



The concept of a star entails a configuration of matter which remains localized. One therefore 
looks for a localized and time-independent matter configuration such that the gravitational field 
is stationary and regular everywhere. As shown in [91J . such a configuration does not exist for a 
real scalar field. But since the stress-energy tensor depends only on the modulus of the scalar field 
and its gradients, one can relax the assumption of time-independence of the scalar while retaining 
a time-independent gravitational field. The key is to assume a harmonic ansatz for the scalar field 



(r,i) =M*Y' 



(33) 



where 4>o is a rea l scalar which is the profile of the star and w is a real constant denoting the 
angular frequency of the phase of the field in the complex plane. 

We consider spherically symmetric, equilibrium configurations corresponding to minimal energy 
solutions while requiring the space-time to be static. In Schwarzschild-like coordinates, the general, 
spherically symmetric, static metric can be written as 



ds 2 



-a (rf dt 2 + a (r) 2 dr 2 + r 2 dfl 2 , 



(34) 



in terms of two real metric functions, a and a. The coordinate r is an areal radius such that 
spheres of constant r have surface area 4irr 2 . For this reason, these coordinates are often called 
polar-areal coordinates. 

The equilibrium equations are obtained by substituting the metric of Eq. (I34[) and the harmonic 
ansatz of Eq. (|33l) into the spherically symmetric EKG system of Eqs. (f2"Tl- l3"2"j) with p = 0, b = 1, 
resulting in three first order partial differential equations (PDEs) 



d r a = -^(« 2 -l) 
2r 



+ 47rGra 2 r 
4irGraa 2 S r r 



d r <$> = - [l + a 2 + 4nGr 2 a 2 (S r r 



r 



dV 
W 2 



a 2 4> . 



(35) 
(36) 

(37) 



Notice that these equations hold for any stress-energy contributions and for a generic type of self- 
potentials F(|</>| 2 ). In order to close the system of Eqs. (I351- I37)) . we still have to prescribe this 
potential. The simplest case admitting localized solutions is the free field case of Eq. (fT2"j) for which 
the potential describes a field with mass m and for which the equations can be written as 



d r a 



2+2 



a 2 4>\ 



(38) 
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o a -l 



a; 



a- 



{1 



^AnGr 

r 

a 2 - 47rGr 2 a 2 m^ 



2 2i2 

m | a <pQ 



$ 2 



} £_(i 



0o a 



(39) 
(40) 



In order to obtain a physical solution of this system, we have to impose the following boundary 
conditions, 



0o (0) 

*(0) 
o(0) 
lim O (r) 



0. 

1, 
0. 



lim a (r) = lim 



a(r) ' 



(41) 
(42) 
(43) 
(44) 

(45) 



which guarantee regularity at the origin and asymptotic flatness. For a given central value of the 
field {0 C }, we need only to adjust the eigenvalue {w} to find a solution which matches the asymp- 
totic behavior of Eqs. (|44H 45 p . This system can be solved as a shooting problem by integrating 
from r = towards the outer boundary r = r ou t ■ Eq. (|39[) is linear and homogeneous in a and one 
is therefore able to rescale the field consistent with Eq. (|45|) . We can get rid of the constants in 
the equations by re-scaling the variables in the following manner 



0o = v47rG0o , r — mr 



t = U) t 



{m/uj)a. 



(46) 



Notice that the form of the metric in Eq. resembles Schwarzschild allowing the association 
= (1 — 2M/r) _1 , where M is the ADM mass of the spacetime. This allow us to define a more 



general mass aspect function 



M(r, t) 



1 



a 2 (r,t) 



which measures the total mass contained in a coordinate sphere of radius r at time t. 
In isotropic coordinates, the spherically symmetric metric can be written as 

ds 2 = -a {Rf dt 2 + if, (i?) 4 (dR 2 + R 2 dQ 2 ) , 



(47) 



(48) 



where tp is the conformal factor. A change of the radial coordinate R = R(r) can transform the 
solution obtained in Schwarzschild coordinates into isotropic ones, in particular 



R(r B 



1 + Va 



dR 
dr 



R 
a — 

r 



(49) 
(50) 



where the first condition is the initial value to integrate the second equation backwards, obtained 
by imposing that far away from the boson star the spacetime resembles Schwarzschild solution. 

As above, boson stars are spherically symmetric solutions of the Eqs. ([381 - 140 [) with asymptotic 
behaviour given by Eqs. (|41H 45 p . For a given value of the central amplitude of the scalar field 
0o (r = 0) = 0c, there exist configurations with some effective radius and a given mass satisfying 
the previous conditions for a different set of n discrete eigenvalues u/™). As n increases, one obtains 
solutions with an increasing number of nodes in 0o . The configuration without nodes is the ground 
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state, while all those with any nodes are excited states. As the number of nodes increases, the 
distribution of the mass as a function of the radius becomes more homogeneous. 

As the amplitude <j) c increases, the stable configuration has a larger mass while its effective 
radius decreases. This trend indicates that the compactness of the boson star increases. However, 
at some point the mass instead decreases with increasing central amplitude. Similar to models 
of neutron stars (see Section 4 of [60] ). this turnaround implies a maximum allowed mass for a 
boson star in the ground state, which numerically was found to be M max = 0.633M| lanck /m. The 
existence of a maximum mass for boson stars is a relativistic effect, which is not present in the 
Newtonian limit, while the maximum of baryonic stars is an intrinsic property. 

Solutions for a few representative boson stars in the ground state are shown in Figure [2] in 
isotropic coordinates. The boson stars becomes more compact for higher values of 4> c , implying 
narrower profiles for the scalar field, larger conformal factors, and smaller lapse functions, as the 
total mass increases. 




G 20 3E 10 20 

r r 



Figure 2: Profiles characterizing static, spherically symmetric boson stars with a few different 
values of the central scalar field (top left). From I143J . 



11 



3 Varieties of Boson Stars 



Quite a number of different flavours of boson stars are present in the literature. They can have 
charge, a fermionic component, or rotation. They can be constructed with various potentials for 
the scalar field. The form of gravity which holds them together can even be modified to, say, 
Newtonian gravity or even no gravity at all (Q-balls). To a certain extent, such modifications are 
akin to varying the equation of state of a normal, fermionic star. Here we briefly review some of 
these variations, paying particular attention to recent work. 

3.1 Self-interaction potentials 

Originally, boson stars were constructed with a free-field potential without any kind of self- 
interaction, obtaining a maximum mass with a dependence M s» Mp lanck /m. This mass, for 
typical masses of bosonic particle candidates, is much smaller than the Chandrasekhar mass 
Mch ~ ^pianck/ 771,2 obtained for fermionic stars, and so they were known as mini-boson stars. 
In order to extend this limit and reach astrophysical masses comparable to the Chandrasekhar 
mass, the potential was generalized to include a self-interaction term that provided an extra pres- 
sure against gravitational collapse. 

Although the first expansion to nonlinear potentials was considered in |162j including fourth 
and sixth power |<^|-terms, a deeper analysis was performed later considering a potential with only 
the quartic term [58] 

v(\cj)\ 2 ) =m 2 H 2 +^H 4 , (51) 

with A a dimensionless coupling constant. Written in terms of a general potential, the EKG 
equations remain the same. The families of gravitational equilibrium can be parametrized by the 
single dimensionless quantity A = A/ (47rGm 2 ). The potential of Eq. (f5"Tj) results in a maximum 
boson star mass that now scales as 

M max « 0.22A 1 / 2 Af Planck /m = (O.lGeV 2 ) M A 1/2 /m 2 (52) 

which is comparable to the Chandrasekhar mass for fermions with mass m/A 1 / 4 [55]. This self- 
interaction therefore allows much larger masses than the mini-boson stars as long as A 3> 1, 
an inequality that may be satisfied even when A <C 1 for reasonable scalar boson masses. The 
maximum mass as a function of the central value of the scalar field is shown in Figure [3] for 
different values of A. The compactness of the most massive stable stars was studied in jS], finding 
an upper bound M/R < 0.16 for A 3> 1. Figure [4] displays this compactness as a function of A 
along with the compactness of a Schwarzschild BH and non-spinning neutron star for comparison. 

Many subsequent papers further analyze the EKG solutions with polynomial, or even more 
general non-polynomial, potentials. One work in particular [196 studied the properties of the 
galactic dark matter halos modeled with these boson stars. They found that a necessary condition 
to obtain stable, compact solutions with an exponential decrease of the scalar field, the series 
expansion of these potentials must contain the usual mass term m 2 \(j)\ 2 . 

More exotic ideas similarly try to include a pressure to increase the mass of BSs. Ref. [5] 
considers a form of repulsive self-interaction mediated by vector mesons within the mean-field 
approximation. However, the authors leave the solution of the fully nonlinear system of the Klein- 
Gordon and Proca equations to future work. Ref. [19] models stars made from the condensation of 
axions, using the semi-relativistic approach with two different potentials. Mathematically this ap- 
proach involves averages such that the equations are equivalent to assuming the axion is constituted 
by a complex scalar field with harmonic time dependence. 

Other generalizations of the potential allow for the presence of nontopological soliton solutions 
even in the absence of gravity, with characteristics quite different than those of the mini-boson 
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Figure 3: (Left) The mass of the boson star as a function of the central value of the scalar field in 
adimensional units a c — \J \nG<\> c . (Right) Maximum mass as a function of A (squares) and the 
asymptotic A — > oo relation of Eq. (|52|) (solid curve). From [58] , 
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Figure 4: The compactness of a stable boson star (black solid line) as a function of the adimensional 
self-interaction parameter A = A/ (47rGm 2 ) . The compactness is shown for the most massive stable 
star (the most compact BS is unstable). This compactness asymptotes for A — > oo to the value 
indicated by the red, dashed line. Also shown for comparison is the compactness of a Schwarzschild 
BH (green dot-dashed line), and the maximum compactness of a non-spinning neutron star (blue 
dotted line). From [9]. 
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stars. In order to obtain these solutions the potential must satisfy two conditions. First, it must 
be a function of \<ft\ 2 to preserve the global U(l) invariance. Second, the potential should have 
an attractive term, bounded from below and positive for \<f>\ — > oo. These conditions imply a 
potential of at least sixth order, a condition that is satisfied by the typical degenerate vacuum 
form [T491I9T1I92] 

v{\cft\ 2 )= m 2 \<ft\ 2 (l-^f) , (53) 

for which the potential has two degenerate minima at ±0o- The case \<j>\ = corresponds to the 
true vacuum state, while \(ft\ = (fto represents the degenerate vacuum state. 

The resulting soliton solution can be split in three different regions. When gravity is negligible, 
the interior solution satisfies tj> ~ </>o , followed by a shell of width 1/m over which <f> changes from 
4>o to zero, and an exterior that is essentially vacuum. This potential leads to a different scaling 
of the mass and the radius than that of the ground state [151] 

M max «M^ laack /(m^) , i? max «M|i anck /(m^). (54) 

There is another type of non-topological soliton star, called Q-stars [157] , which also admits 
soliton solutions in the absence of gravity (i.e., Q-balls [5711151] ). The potential, besides being also 
a function of \(f>\ 2 , must satisfy the following conditions: it must behave like ~ (ft 2 near (ft = 0, it has 
to be bounded < (j) 2 in an intermediate region and must be larger > (ft 2 for \(ft\ — > oo. The Q-stars 
also have three regions; an interior solution of radius R s=s Mpi anc k/</>§, (i.e., (fto « m is the free 
particle inverse Compton wavelength) a very thin surface region of thickness I /(fto, and finally the 
exterior solution without matter, which reduces to Schwarzschild in spherical symmetry. The mass 
of these Q-stars scales now as Mp lanck /(pQ. The stability of these Q-stars has been studied recently 
using catastrophe theory, such as [210| 1137] . Rotating, axisymmetric Q-balls were constructed in 
|135| 1136] . Related, rotating solutions in 2 + 1 with the signum-Gordon equation instead of the 
KG equation are found in [TTj . 

Other interesting works have studied the formation of Q-balls by the Affleck-Dine mecha- 
nism |127j . their dynamics in one, two and three spatial dimensions [23], and their viability as a 
self-interacting dark matter candidate [141] . 

Ref. [30] considers a chemical potential to construct BSs, arguing that the effect of the chemical 
potential is to reduce the parameter space of stable solutions. Related work modifies the kinetic 
term of the action instead of the potential. Ref. [5] studies the resulting BSs for a class of K 
field theories, finding solutions of two types: (i) compact balls possessing a naked singularity at 
their center and (ii) compact shells with a singular inner boundary which resemble black holes. 
Ref. 14 considers coherent states of a scalar field instead of a BS within k-essence in the context of 
explaining dark matter. Ref. [73] modifies the kinetic term with just a minus sign to convert the 
scalar field to a phantom field. Although, a regular real scalar field has no spherically symmetric, 
local static solutions, they find such solutions with a real phantom scalar field. 

3.2 Newtonian boson stars 

The Newtonian limit of the Einstein-Klein-Gordon Eqs. (I91 ITT]) can be derived by assuming that 
the spacetime metric in the weak field approximation can be written as 

g ao = -(l + 2V) , 9ii = l + 2V, gij =0 for i^j, (55) 

where V is the Newtonian gravitational potential. In this limit, the Einstein equations reduce to 
the Poisson equation 

V 2 V = 47tGT 00 = 4irGm 2 (ft4>. (56) 
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Conversely, by assuming that 

4>(x,t) = ^{x,t)e im \ (57) 
in addition to the weak limit of Eq. (|55[) , the Klein-Gordon equation reduces to 

id t ® = -V 2 $ + my$ (58) 

2 m 

which is just the Schrodinger equation with h = 1. Therefore, the EKG system is reduced in the 
Newtonian limit to the Schrodinger-Poisson (SP) system [100 . 

The initial data is obtained by solving an eigenvalue problem very similar to the one for boson 
stars, with similar assumptions and boundary conditions. The solutions also share similar features 
and display a similar behaviour. A nice property of the Newtonian limit is that all the solutions 
can be obtained by rescaling from one known solution jlOOj . 

02=01 (§9 ' w2=cji GI) ' r2=ri GI)' (59) 

where N = m J dx 3 ^ is the Newtonian number of particles. 

The possibility of including self-interaction terms in the potential was considered in [108] , study- 
ing also the gravitational cooling (i.e., the relaxation and virialization through the emission of scalar 
field bursts) of spherical perturbations. Non-spherical perturbations were further studied in [28] , 
showing that the final state is a spherically symmetric configuration. Single Newtonian boson stars 
were studied in [100] , either when they are boosted with/without an external central potential. Ro- 
tating stars were first successfully constructed in |204j within the Newtonian approach. Numerical 
evolutions of binary boson stars in Newtonian gravity are discussed in Section 14.21 

Recent work by Chavanis with Newtonian gravity solves the Gross-Pitaevskii equation, a vari- 
ant of Eq. (|58[) which involves a pseudo-potential for a Bose-Einstein condensate, to model either 
dark matter or compact alternatives to neutron stars [46j [47] 08] . 

Much recent work considers boson stars from a quantum perspective as a Bose-Einstein con- 
densate involving some number, N, of scalar fields. Ref. [161] studies the collapse of boson stars 
mathematically in the mean field limit in which N —> oo. Ref. |132] argues for the existence of 
bosonic atoms instead of stars. Ref. [17] uses numerical methods to study the mean field dynamics 
of BSs. 



3.3 Charged boson stars 

Charged boson stars result from the coupling of the boson field to the electromagnetic field [125] , 
The coupling between gravity and a complex scalar field with a U(l) charge arises by considering 
the action of Eq. [7] with the following matter Lagrangian density 



g ab (Va4> + ieA a 4>)(V b <j)-ieA b <b) + v(\cP\ 2 ^] -^F ab F ab , (60) 



where e is the gauge coupling constant. The Maxwell tensor F ab can be decomposed in terms of 
the vector potential A a 

F ab = V a A b - V b A a . (61) 

The system of equations obtained by performing the variations on the action forms the Einstein- 
Maxwell-Klein-Gordon system, which contains the evolution equations for the complex scalar field 
<j>, the vector potential A a and the spacetime metric g ab [180j . 

Because a charged BS may be relevant for a variety of scenarios, we detail the resulting equa- 
tions. For example, cosmic strings are also constructed from a charged, complex scalar field and 
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obeys these same equations. It is only when we choose the harmonic time dependence of the scalar 
field that we distinguish from the harmonic azimuth of the cosmic string |219j . The evolution 
equations for the scalar field and for the Maxwell tensor are 



g ab V a V b (t) - 2 i eA 



V a (f>-e 2 <f>A a A a -ie 



<t>V a A 



dV 



(62) 



y aF ab = _jb = ie (0 Vfc _ Vfe 0) + 2 e 2 (j)A b . 



(63) 



Notice that the vector potential is not unique; we can still add any curl-free components without 
changing the Maxwell equations. The gauge freedom can be fixed by choosing, for instance, the 
Lorentz gauge V ' a A a — 0. Within this choice, which sets the first time derivative of the time 
component Aq 1 the Maxwell equations reduce to a set of wave equations in a curved background 
with a non-linear current. This gauge choice resembles the harmonic gauge condition, which casts 
the Einstein equations as a system of non-linear, wave equations |220j . 

Either from Noether's theorem or by taking an additional covariant derivative of Eq. (|63p . one 
obtains that the electric current J a follows a conservation law. The spatial integral of the time 
component of this current, which can be identified with the total charge Q, is conserved. This 
charge is proportional to the number of particles, Q = e N. The mass M and the total charge 
Q can be calculated by associating the asymptotic behavior of the metric with that of Reissner- 
Nordstrom metric, 



which is the unique solution at large distances for a scalar field with compact support. 

We look for a time independent metric by first assuming a harmonically varying scalar field 
as in Eq. (|33p . We work in spherical coordinates and assume spherical symmetry. With a proper 
gauge choice, the vector potential takes a particularly simple form with only a single, non-trivial 
component A a — (Ao(r), 0, 0, 0). This choice implies an everywhere vanishing magnetic field so that 
the electromagnetic field is purely electric. The boundary conditions for the vector potential are 
obtained by requiring the electric field to vanish at the origin because of regularity, d r A (r = 0) = 0. 
Because the electromagnetic field depends only on derivatives of the potential, we can use this 
freedom to set A (oo) = [125] . 

With these conditions, it is possible to find numerical solutions in equilibrium as described in 
Ref. [125] . Solutions are found for e 2 = e 2 M"pi anck /(8 n m 2 ) < 1/2. For e 2 > 1/2 the repulsive 
Coulomb force is bigger than the gravitational attraction and no solutions are found. This bound 
on the BS charge in terms of its mass ensures that one cannot construct an overcharged BS, 
in analogy to the overcharged monopoles of Ref. |156j . An overcharged monopole is one with 
more charge than mass and is therefore susceptible to gravitational collapse by accreting sufficient 
(neutral) mass. However, because its charge is higher than its mass, such collapse might lead to an 
extremal Reissner-Nordstrom BH, but BSs do not appear to allow for this possibility. Interestingly 
Ref. |191] finds that if one removes gravity, the obtained Q-balls may have no limit on their charge. 

The mass and the number of particles are plotted as a function of <j> c for different values of 
e in Figure [SJ Trivially, for e = the mini-boson stars of Section (|2.4p are recovered. Excited 
solutions with nodes are qualitatively similar (1 25] . The stability of these objects has been studied 
in [121] , showing that the equilibrium configurations with a mass larger than the critical mass are 
dynamically unstable, similar to uncharged BSs. 

Recent work with charged BSs includes the publication of Maple [1] routines to study boson 
nebulae charge [Ml EHl ESS] 

Other work generalizes the Q-balls and Q-shells found with a certain potential which leads to 
the signum-Gordon equation for the scalar field |133l 1134] , In particular, shell solutions can be 




for r — > oo 



(64) 
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Figure 5: The mass (solid) and the number of particles (dashed) versus central scalar value for 
charged boson stars with four values of e as defined in Section 13.31 The mostly- vertical lines 
crossing the four plots indicate the solution for each case with the maximum mass (solid) and 
maximum particle number (dashed). From |125j . 



found with a black hole in its interior which has implications for black hole scalar hair (for a review 
of black hole uniqueness see |115j ). 

One can also consider Q-balls coupled to an electromagnetic field, a regime appropriate for 
particle physics. Within such a context, Ref. [77] studies the chiral magnetic effect arising from a 
Q-ball. Other work in Ref. [37] studies charged, spinning Q-balls. 



3.4 Oscillatons 

As mentioned earlier, it is not possible to find time- independent, spacetime solutions for a real 
scalar field. However, there are non-singular, time-dependent near-equilibrium configurations of 
self-gravitating real scalar fields, which are known as oscillatons |198j . These solutions are similar 
to boson stars, with the exception that the spacetime must also have a time dependence in order 
to avoid singularities. 

In this case, the system is still described by the EKG Eqs. (|27H32j) . with the the additional 
simplification that the scalar field is strictly real, (f> = <p. In order to find equilibrium configurations, 
one expands both metric components {A(r, t) = a 2 ,C(r,t) = (a/a) 2 } and the scalar field <f>(r,i) as 
a truncated Fourier series 

im»x 

<f>(r,t) = 5^_i(r)coB([2j-l]wi) , (65) 
i=i 

jmax Jmax 

A(r,t) = Y,A 2j (r)cos(2jut) , C{r, t) = £ C 2i (r) cos(2jwt) , (66) 

3=0 3=0 

where w is the fundamental frequency and j max is the mode at which the Fourier series are trun- 
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cated. As noted in Ref. [217) 16]. the scalar field consists only of odd components while the metric 
terms consist only of even ones. Solutions are obtained by substituting the expansions of Eq. (|65|) 
into the spherically symmetric Eqs. (|27M32[) . By matching terms of the same frequency, the system 
of equations reduce to a set of coupled ODEs. The boundary conditions are determined by requiring 
regularity at the origin and that the fields become asymptotically flat at large radius. These form 
an eigenvalue problem for the coefficients {^j-ifV = 0),A 2 j(r = 0),C2j(r = 0)} corresponding 
to a given central value <j>i(r = 0). As pointed out in |217j . the frequency u> is determined by the 
coefficient Co(oo) and is therefore called an output value. Although the equations are non-linear, 
the Fourier series converges rapidly, and so a small value of j ma x usually suffices. 

A careful analysis of the high frequency components of this construction reveals difficulties in 
avoiding infinite total energy while maintaining the asymptotically flat boundary condition [173] , 
Therefore, the truncated solutions constructed above are not exactly time periodic. Indeed, very 
accurate numerical work has shown that the oscillatons radiate scalar field on extremely long time 
scales while their frequency increases [86] [99]. This work finds a mass loss rate of just one part 
in 10 12 per oscillation period, much too small for most numerical simulations to observe. The 
solutions are therefore only near-equilibrium solutions and can be extremely long-lived. 
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Figure 6: (Top) Total mass (in units of ^pi anc k/ m ) an d fundamental frequency of an oscillaton 
as a function of the central value of the scalar field (j>i(r — 0). The maximum mass is M max = 
0.607M| lanck /m. (Bottom) Plot of the total mass versus the radius at which g rr achieves its 
maximum. From [BJ. 

Although the geometry is oscillatory in nature, these oscillatons behave similar to BSs. In 
particular, they similarly transition from long-lived solutions to a dynamically unstable branch 
separated at the maximum mass M max = 0.607 Mp lanck /m. Figure [SJ displays the total mass 
curve which shows the mass as a function of central value. Compact solutions can be found 
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in the Newtonian framework when the weak field limit is performed appropriately, reducing to 
the so-called Newtonian oscillations [217] . The dynamics produced by perturbations are also 
qualitatively similar, including gravitational cooling, migration to more dilute stars, and collapse 
to black holes [6]. More recently, these studies have been extended by considering the evolution 
in 3D of excited states [14] and by including a quartic self- interaction potential [218] . In [131] . 
a variational approach is used to construct oscillatons in a reduced system similar to that of the 
sine-Gordon breather solution. 

Closely related, are oscillons that exist in flatspace and that were first mentioned as "pulsons" 
in 1977 [34] , There is an extensive literature on such solutions, many of which appear in [82] , A 
series of papers establishes that oscillons similarly radiate on very long time scales (82] [83] [84] [85] . 
An interesting numerical approach to evolving oscillons adopts coordinates that blueshift and 
damp outgoing radiation of the massive scalar field |117[ 1119] . A detailed look at the long term 
dynamics of these solutions suggests the existence of a fractal boundary in parameter space between 
oscillatons that lead to expansion of a true- vacuum bubble and those that disperse |118] . 



3.5 Rotating boson stars 

Boson stars with rotation were not explored until the mid-1990s because of the lack of a strong 
astrophysical motivation and the technical problems with the regularization along the axis of sym- 
metry. The first equilibrium solutions of rotating boson stars were obtained within the Newtonian 
gravity [204] . In order to generate axisymmetric time-independent solutions with angular momen- 
tum, one is naturally lead to the ansatz 

cj)(r,t) = Mr,dy {ut+k *\ (67) 

where </>o( r , 0) is a real scalar representing the profile of the star, w is a real constant denoting the 
angular frequency of the field and k must be an integer so that the field cj> is not multivalued in 
the azimuthal coordinate ip. This integer is commonly known as the rotational quantum number. 

General relativistic rotating boson stars were later found |194| 1223] with the same ansatz of 
Eq. (]67[) . To obtain stationary axially symmetric solutions, two symmetries were imposed on the 
spacetime described by two commuting Killing vector fields £ = dt and rj = d v in a system of 
adapted (cylindrical) coordinates {t,r,0,tp}. In these coordinates, the metric is independent of t 
and ip and can be expressed in isotropic coordinates in the Lewis-Papapetrou form 



ds 2 = -fdt 2 + - 



g (dr 2 + r 2 dO 2 ) + r 2 sin 2 9 [dip- -dt 



where f,l,g and fl are metric functions depending only on r and 9. This means that we have 
to solve five coupled PDEs, four for the metric and one for the Klein-Gordon equation; these 
equations determine an elliptic quadratic eigenvalue problem in two spatial dimensions. Near the 
axis, the scalar field behaves as 

lim Mr, 0) = r k h k {9) + 0{r r+2 ) (69) 

r— >0 

so that for k > the field vanishes near the axis. Note that hk is some arbitrary function different 
for different values of k but no sum over k is implied in Eq. (I69[) . This implies that the rotating 
star solutions have toroidal level surfaces instead of spheroidal ones as in the spherically symmetric 
case k = 0. In this case the metric coefficients are simplified, namely g = 1, CI = and / = f(r), 
l = l(r). 

The entire family of solutions for k = 1 and part of k = 2 was computed using the self- 
consistent field method [223], obtaining a maximum mass M max = 1.31Mp lanck /m. Both families 
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were completely computed in [143] using faster multigrid methods, although there were significant 
discrepancies in the maximum mass which indicates a problem with the regularity condition on the 
z-axis. The mass M and angular momentum J for stationary asymptotically flat spacetimes can 
be obtained from their respective Komar expressions. They can be read off from the asymptotic 
expansion of the metric functions / and 

+ n-^ + o(£). <™> 

Alternatively, using the Tolman expressions for the angular momentum and the Noether charge 
relation in Eq. ll5[ one obtains an important quantization relation for the angular momentum [223 

J = kN, (71) 

for integer values of k. This remarkable quantization condition for this classical solution also plays 
a role in the work of [70] discussed in Section 16.31 Figure [7] shows the scalar field for two different 
rotating BSs. 

Recently, their stability properties were found to be similar to nonrotating stars [137] . Rotating 
boson stars have been shown to develop a strong ergoregion instability when rapidly spinning on 
short characteristic timescales (i.e., 0.1 seconds- 1 week for objects with mass M = 1 — 10 6 M Q 
and angular momentum J > OAGM 2 ), indicating that very compact objects with large rotation 
are probably black holes [4"5] . 

More discussion concerning the numerical methods and limitations of some of these approaches 
can be found in Lai's PhD thesis [143] . 



3.6 Fermionic-bosonic stars 

The possibility of compact stellar objects made with a mixture of bosonic and fermionic matter 
was studied in Refs. (1 131 1114J . In the simplest case, the boson component interacts with the 
fermionic component only via the gravitational field, although different couplings were suggested 
in [114] and have been further explored in [67 ]ll81j . Such a simple interaction is, at the very least, 
consistent with models of a bosonic dark matter coupling only gravitationally with visible matter, 
and the idea that such a bosonic component would become gravitationally bound within fermionic 
stars is arguably a natural expectation. 

One can consider a perfect fluid as the fermionic component such that the stress-energy tensor 
takes the standard form 



T, 



fluid 



(fi + p)u a u b +pg ab 



(72) 



where y, is the energy density, p is the pressure of the fluid and u a its four-velocity. Such a fluid 
requires an equation of state to close the system of equations (see Ref. [57] for more about fluids 
in relativity). In most work with fermionic-boson stars, the fluid is described by a degenerate, 
relativistic Fermi gas, so that the pressure is given by the parametric equation of state of Chan- 
drasekhar 



|U = K (sinht — t) 



K 



sinh t 







G) 


+ 3t 





(73) 



where K = m^/ (32ir 2 ) and TTt n the mass of the fermion. The parameter t is given by 



t(r) = 4 log 




(74) 
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Figure 7: The scalar field in cylindrical coordinates <p(p, z) for two rotating boson star solutions: 
(top) k = 1 and (bottom) k — 2. The two solutions have roughly comparable amplitudes in scalar 
field. Note the toroidal shape. From [143j . 
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where p is the maximum value of the momentum in the Fermi distribution at radius r. 

The perfect fluid obeys relativistic versions of the Euler equations, which account for the con- 
servation of the fluid energy and momentum, plus the conservation of the baryonic number (i.e., 
mass conservation). The complex scalar field representing the bosonic component is once again 
described by the Klein-Gordon equation. The spacetime is computed through the Einstein equa- 
tions with a stress-energy tensor which is a combination of the complex scalar field and the perfect 
fluid 



Tab = T, 



"i fluid 



ab "t" ab 



(75) 



After imposing the harmonic time dependence of Eq. (|33[) on the complex scalar field, assuming 
a static metric as in Eq. (]34l) and the static fluid u a = 0, one obtains the equations describing 
equilibrium fermion-boson configurations 
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These equations can be written in adimensional form by rescaling the variables by introducing 
the following quantities 

x = mr , a(x) = v / 47rG0(O, r) , ee cu/m 2 , /j ee (AirG/m 2 )^ , p ee (47tG/to 2 )p. (76) 

By varying the central value of the fermion density fi(r = 0) and the scalar field <f>(r — 0), one 
finds stars dominated by either bosons or fermions, with a continuous spectrum in between. 

It was shown that the stability arguments made with boson stars can also be applied to these 
mixed objects [123J . The existence of slowly rotating fermion-boson stars was shown in [66] . 
although no solutions were found in previous attempts [138] , Also see [74] for unstable solutions 
consisting of a real scalar field coupled to a perfect fluid with a polytropic equation of state. 



3.7 Multi-state boson stars 

It turns out that excited BSs, as dark matter halo candidates, provide for flatter, and hence more 
realistic, galactic rotation curves than ground state BSs. The problem is that they are generally 
unstable to decay to their ground state. Combining excited states with the ground states in what 
are called multi-state BSs is one way around this. 

Although bosons in the same state are indistinguishable, it is possible to construct non-trivial 
configurations with bosons in different excited states. A system of bosons in P different states that 
only interact with each other gravitationally can be described by the following Lagrangian density 



C 



1 

16ttG 



R 



^ 2 

n=l 



(77) 
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where is the particular complex scalar field representing the bosons in the n-state with n — 1 
nodes. The equations of motion are very similar to the standard ones described in Section |2~21 with 
two peculiarities: (i) there are n independent KG equations (i.e., one for each state) and (ii) the 
stress-energy tensor is now the sum of contributions from each mode. Equilibrium configurations 
for this system were found in [26] . 

In the simplest case of a multi-state boson, one has the ground state and the first excited state. 
Such configurations are stable if the number of particles in the ground state is larger than the 
number of particles in the excited state H] 

N {1) >N {2) . (78) 

This result can be understood as the ground state deepens the gravitational potential of the excited 
state, and thereby stabilizing it. Unstable configurations migrate to a stable one via a flip-flop of 
the modes; the excited state decays, while the ground jumps to the first exited state, so that the 
condition (|78l) is satisfied. An example of this behaviour can be observed in Figure [SJ 




Figure 8: (Top) The maximum of the central value of each of the two scalar fields constituting the 
multi-state BS for the fraction r\ = 3 where rj = /N^ defines the relative "amount" of each 
state. (Bottom) The frequencies associated with each of the two states of the multi-state BS. At 
t = 2000, there is a flip in which the excited state (black solid) decays and the scalar field in the 
ground state (red dashed) becomes excited. Discussed in Section |3~T1 From [26 . 

Similar results were found in the Newtonian limit 216 a , however with a slightly higher stability 
limit > 1.13 N^ 2 \ This work stresses that combining several excited states makes it possible to 
obtain flatter rotation curves than only with ground state, producing better models for galactic dark 
matter halos (see also discussion of boson stars as an explanation of dark matter in Section HT3")) . 
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Ref. |112] considers two scalar fields describing boson stars which are phase shifted in time with 
respect to each other, studying the dynamics numerically. In particular, one can consider multiple 
scalar fields with an explicit interaction (beyond just gravity) between them, say V (l*^ 1 ^ |<?^ 2 ^|)- 
Refs. [38] [39] construct such solutions, considering the individual particle-like configurations for 
each complex field as interacting with each other. 

3.8 Alternative theories of gravity 

Instead of modifying the scalar field potential, one can instead consider alternative theories of 
gravity. Constraints on such theories are already significant given the great success of general 
relativity |222j . However, the fast advance of electromagnetic observations and the anticipated 
gravitational wave observations promise much more in this area, in particular in the context of 
compact objects that probe strong-field gravity. 

An ambitious effort is begun in Ref. [177] which studies a very general gravitational Lagrangian 
("extended scalar-tensor theories") with both fluid stars and boson stars. The goal is for observa- 
tions of compact stars to constrain such theories of gravity. 

It has been found that scalar tensor theories allow for spontaneous scalarization in which the 
scalar component of the gravity theory transitions to a non-trivial configuration analogously to 
ferromagnetism with neutron stars [63] . Such scalarization is also found to occur in the context of 
boson star evolution [7]. 

Boson stars also occur within conformal gravity and with scalar-tensor extensions to it [411 142] . 

One motivation for alternative theories is to explain the apparent existence of dark matter 
without resorting to some unknown dark matter component. Perhaps the most well known of 
these is MOND (modified Newtonian dynamics) in which gravity is modified only at large dis- 
tances [165| 1166] (for a review see [78]). Boson stars are studied within TeVeS (Tensor- Vector- 
Scalar), a relativistic generalization of MOND [59]. In particular, their evolutions of boson stars 
develop caustic singularities, and the authors propose modifications of the theory to avoid such 
problems. 

3.9 Gauged boson stars 

In 1988, Bartnik and McKinnon published quite unexpected results showing the existence of 
particle-like solutions within SU{2) Yang-Mills coupled to gravity [21]. These solutions, although 
unstable, were unexpected because no particle-like solutions are found in either the Yang-Mills 
or gravity sectors in isolation. Recall also that no particle-like solutions were found with gravity 
coupled to electromagnetism in early efforts to find Wheeler's geon (however, see Section 16.31 for 
discussion of Ref. 71 which finds geons within AdS). 

Bartnik and McKinnon generalize from the Abelian U{\) gauge group to the non-Abelian 
SU{2) group and thereby find these unexpected particle-like solutions. One can consider, as does 
Ref. [195 j (see Section Hp), these globally regular solutions (and their generalizations to SU(n) for 
n > 2) as gauged boson stars even though these contain no scalar field. One can instead explicitly 
include a scalar field doublet coupled to the Yang-Mills gauge field [40] as perhaps a more direct 
generalization of the (U(l)) charged boson stars discussed in Section [3T3T 

Ref. [75] studies BSs formed from a gauge condensate of an SU(3) gauge field, and Ref. [42] 
extends the Bartnik-McKinnon solutions to conformal gravity with a Higgs field 42 . 
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4 Dynamics of Boson Stars 



In this section, the formation, stability and dynamical evolution of boson stars are discussed. One 
approach to the question of stability considers small perturbations around an equilibrium config- 
uration, so that the system remains in the linearized regime. Growing modes indicate instability. 
However, a solution can be linearly stable and yet have a nonlinear instability. One example 
is Minkowski space which, under small perturbations, relaxes back to flat, but, for sufficiently 
large perturbations, leads to black hole formation, decidedly not Minkowski. To study nonlinear 
stability, other methods are needed. In particular, full numerical evolutions of the Einstein-Klein- 
Gordon (EKG) equations are quite useful for understanding the dynamics of boson stars. 

4.1 Gravitational stability 

A linear stability analysis consists of studying the time evolution of infinitesimal perturbations 
about an equilibrium configuration, usually with the additional constraint that the total number 
of particles must be conserved. In the case of spherically symmetric, fcrmionic stars described by 
perfect fluid, it is possible to find an eigenvalue equation for the perturbations that determines 
the normal modes and frequencies of the radial oscillations (see, for example, Ref. [55]). Stability 
theorems also allow for a direct characterization of the stability branches of the equilibrium solu- 
tions [S3] [HI] ■ Analogously, one can write a similar eigenvalue equation for boson stars and show 
the validity of similar stability theorems. In addition to these methods, the stability of boson stars 
has also been studied using two other, independent methods: by applying catastrophe theory and 
by solving numerically the time dependent Einstein-Klein-Gordon equations. All these methods 
agree with the results obtained in the linear stability analysis. 

4.1.1 Linear stability analysis 

Assume that a spherically symmetric boson star in an equilibrium configuration is perturbed only 
in the radial direction. The equations governing these small radial perturbations are obtained by 
linearizing the system of equations in the standard way; expand the metric and the scalar field 
functions to first order in the perturbation and neglect higher order terms in the equations [9511120] , 
Considering the collection of fields for the system /j, one expands them in terms of the background 
solution fi and perturbation as 

fi(r,t) = °f i (r) + 1 f i (r)e i ° t (79) 

which assumes harmonic time dependence for the perturbation. Substitution of this expansion 
into the system of equations then provides a linearized system which reduces to a set of coupled 
equations that determines the spectrum of modes 1 fi and eigenvalues a 2 

Uj 1 fi = a 2 Mij 1 fi (80) 

where Ljj is a differential operator containing partial derivatives and My is a matrix depending 
on the background equilibrium fields Solving this system, known as the pulsation equation, 
produces the spectrum of eigenmodes and their eigenvalues a. 

The stability of the star depends crucially on the sign of the smallest eigenvalue. Because of 
time reversal symmetry, only a 2 enters the equations [150] . and we label the smallest eigenvalue 
<7q. If it is negative, the eigenmode grows exponentially with time and the star is unstable. On 
the other hand, for positive eigenvalues the configuration has no unstable modes and is therefore 
stable. The critical point at which the stability transitions from stable to unstable therefore occurs 
when the smallest eigenvalue vanishes, a a = 0. 
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Equilibrium solutions can be parametrized with a single variable, such as the central value of 
the scalar field <f> c , and so we can write M = M(<j> c ) and N = N(<fi c ). Stability theorems then 
indicate that transitions between stable and unstable configurations occur only at critical points 
of the parameterization (M'(<p c ) = N'(<j> c ) — 0) [6T1 l93l [T09l 1208] . Linear perturbation analysis 
provides a more detailed picture such as the growth rates and the eigenmodes of the perturbations. 

Ref. [96^ carries out such an analysis for perturbations that conserve mass and charge. They 
find the first three perturbative modes and their growth rates, and they identify at which precise 
values of (j) c these modes become unstable. Starting from small values, they find that ground state 
BSs are stable up to the critical point of maximum mass. Further increases in the central value 
subsequently encounter additional unstable modes. This same type of analysis applied to excited 
state BSs showed that the same stability criterion applies for perturbations that conserve the 
total particle number 122 . For more general perturbations that do not conserve particle number, 
excited states are generally unstable to decaying to the ground state. 

A more involved analysis by [150] uses a Hamiltonian formalism to study BS stability. Consid- 
ering first order perturbations that conserve mass and charge (SN = 0), their results agree with 
those of 96, 122 . However, they extend their approach to consider more general perturbations 
which do not conserve the total number of particles (i.e., SN ^ 0). To do so, they must work with 
the second order quantities. They found complex eigenvalues for the excited states that indicate 
that excited state boson stars are unstable. More detail and discussion on the different stability 
analysis can be found in Ref. [124] , 

Catastrophe theory is part of the study of dynamical systems that began in the 1960s and 
studies large changes in systems resulting from small changes to certain important parameters (for 
a physics-oriented review see (205 ) . Its use in the context of boson stars is to evaluate stability, 
and to do so one constructs a series of solutions in terms of a limited and appropriate set of 
parameters. Under certain conditions, such a series generates a curve smooth everywhere except 
for certain points. Within a given smooth expanse between such singular points, the solutions 
share the same stability properties. In other words, bifurcations occur at the singular points so 
that solutions after the singularity gain an additional, unstable mode. Much of the recent work in 
this area confirms the previous conclusions from linear perturbation analysis |2101 1211[ 12121 1213] 
and from earlier work with catastrophe theory [142J . 

Another recent work using catastrophe theory finds that rotating stars share a similar stability 
picture as nonrotating solutions [137 . However, only fast spinning stars are subject to an ergoregion 
instability [45] . 



4.1.2 Non-linear stability of single boson stars 

The dynamical evolution of spherically symmetric perturbations of boson stars has also been stud- 
ied by solving numerically the Einstein-Klein-Gordon equations (Section 12 . 3[) . or its Newtonian 
limit (Section 13. 2p . the Schrodinger-Poisson system. The first such work was Ref. [197] in which 
the stability of the ground state was studied by considering finite perturbations which may change 
the total mass and the particle number (i.e., SN ^ and 5M ^ 0). The results corroborated 
the linear stability analysis in the sense that they found a stable and an unstable branch with a 
transition between them at a critical value, </> cr it, of the central scalar field corresponding to the 
maximal BS mass M max = 0.633Afp lanck /rn. 

The perturbed configurations of the stable branch may oscillate and emit scalar radiation 
maintaining a characteristic frequency u, eventually settling into some other stable state with less 
mass than the original. This characteristic frequency can be approximated in the non-relativistic 
limit as [T97] 

7r mGM 

v = ~ (81) 
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where R is the effective radius of the star and M its total mass. Scalar radiation is the only damping 
mechanism available because spherical symmetry does not allow for gravitational radiation and 
because the Klein-Gordon equation has no viscous or dissipative terms. This process was named 
gravitational cooling, and it is extremely important in the context of formation of compact bosonic 
objects [199] (see below). The behavior of perturbed solutions can be represented on a plot of 
frequency versus effective mass as in Figure [9j Perturbed stars will oscillate with a frequency 
below its corresponding solid line and they radiate scalar field to infinity. As they do so, they 
lose mass by oscillating at constant frequency, moving leftward on the plot until they settle on the 
stable branch of (unperturbed) solutions. 




Figure 9: Oscillation frequencies of various boson stars are plotted against their mass. Also shown 
are the oscillation frequencies of unstable BSs obtained from the fully nonlinear evolution of the 
dynamical system. Unstable BSs are observed maintaining a constant frequency as they approach 
a stable star configuration. From [197] , 

The perturbed unstable configurations will either collapse to a black hole or migrate to a stable 
configuration, depending on the nature of the initial perturbation. If the density of the star is 
increased, it will collapse to a black hole. On the other hand, if it is decreased, the star explodes, 
expanding quickly as it approaches the stable branch. Along with the expansion, energy in the 
form of scalar field is radiated away, leaving a very perturbed stable star, less massive than the 
original unstable one. 

This analysis was extended to boson stars with self-interaction and to excited BSs in Ref. [I6j . 
showing that both branches of the excited states were intrinsically unstable under generic pertur- 
bations that do not preserve M and N. The low density excited stars, with masses close to the 
ground state configurations, will evolve to ground state boson stars when perturbed. The more 
massive configurations form a black hole if the binding energy Eb = M — Nm is negative, through 
a cascade of intermediate states. The kinetic energy of the stars increases as the configuration gets 
closer to Eb — 0, so that for positive binding energies there is an excess of kinetic energy that 
tends to disperse the bosons to infinity. These results are summarized in Figure I1Q[ which shows 
the time scale of the excited star to decay to one of these states. 
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Instability Time Scale of a 1-node Star 
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Figure 10: The instability time scale of an excited boson star (the first excitation) to one of three 
end states: (i) decay to the ground state, (ii) collapse to a black hole, or (hi) dispersal. From [16] . 



More recently, the stability of the ground state was revisited with 3D simulations using a Carte- 
sian grid [104 . The Einstein equations were written in terms of the BSSN formulation [2011 124] , 
which is one of the most commonly used formulations in numerical relativity. Intrinsic numerical 
error from discretization served to perturb the ground state for both stable and unstable stars. It 
was found that unstable stars with negative binding energy would collapse and form a black hole, 
while ones with positive binding energy would suffer an excess of kinetic energy and disperse to 
infinity. 

That these unstable stars would disperse, instead of simply expanding into some less compact 
stable solution, disagrees with the previous results of Ref. [197] . and was subsequently further 
analyzed in [106] in spherical symmetry with an explicit perturbation (i.e., a Gaussian shell of 
particles which increases the mass of the star around 0.1%). The spherically symmetric results 
corroborated the previous 3D calculations, suggesting that the slightly perturbed configurations 
of the unstable branch have three possible endstates: (i) collapse to BH, (ii) migration to a less 
dense stable solution, or (iii) dispersal to infinity, dependent on the sign of the binding energy. 

Closely related is the work of Lai and Choptuik [144] studying BS critical behavior (discussed in 
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Figure 11: Very long evolutions of a perturbed, slightly sub-critical, boson star with differing outer 
boundaries. The central magnitude of the scalar field is shown. At early times (t < 250 and the 
middle frame), the boson star demonstrates near-critical behavior with small-amplitude oscillations 
about an unstable solution. For late times (t > 250), the solution appears converged for the largest 
two outer boundaries and suggests that sub-critical boson stars are not dispersing. Instead, they 
execute large amplitude oscillations about low-density boson stars. From |144j . 

Section |6~T|) . They tune perturbations of boson stars so that dynamically the solution approaches 
some particular unstable solution for some finite time. They then study evolutions which ultimately 
do not collapse to BH, so-called sub-critical solutions, and find that they do not disperse to infinity, 
instead oscillating about some less compact, stable star. They show results with increasingly distant 
outer boundary that suggest that this behavior is not a finite-boundary-related effect (reproduced 
in Figure [TTJ . They use a different form of perturbation than Ref. |106j , and, being only slightly 
subcritical, may be working in a regime with non-positive binding energy. However, it is interesting 
to consider that if indeed there are three distinct end-states, then one might expect critical behavior 
in the transition among the different pairings. Non-spherical perturbations of boson stars have been 
studied numerically in [15] with a 3D code to analyze the emitted gravitational waves. 

Much less is known about rotating BSs which are more difficult to construct and to evolve 
because they are, at most, axisymmetric, not spherically symmetric. However, as mentioned in 
Section |3"31 they appear to have both stable and unstable branches [137] and are subject to an 
ergoregion instability at high rotation rates [IS]. To our knowledge, no one has evolved rotating 
BS initial data. However, as discussed in the next section, simulations of BS binaries [168, 174 
have found rotating boson stars as a result of merger. 

The issue of formation of boson stars has been addressed in [199] by performing numerical 
evolutions of the EKG system with different initial Gaussian distributions describing unbound 
states (i.e., the kinetic energy is larger than the potential energy). Quite independent of the initial 
condition, the scalar field collapses and settles down to a bound state by ejecting some of the 
scalar energy during each bounce. The ejected scalar field carries away excess ever-decreasing 
amounts of kinetic energy, as the system becomes bounded. After a few free-fall times of the initial 
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configuration, the scalar field has settled into a perturbed boson star on the stable branch. This 
process is the already mentioned gravitational cooling, and allows for the formation of compact 
soliton stars (boson stars for complex scalar fields and oscillatons for real scalar fields). Although 
these evolutions assumed spherical symmetry which does not include important processes such 
as fragmentation or the formation of pancakes, they demonstrate the feasibility of the formation 
mechanism; clouds of scalar field will collapse under their own self-gravity while shedding excess 
kinetic energy. The results also confirm the importance of the mass term in the potential. By 
removing the massive term in the simulations, the field collapses, rebounds and completely disperses 
to infinity, and no compact object forms. The evolution of the scalar field with and without the 
massive term is displayed in Figure 1121 
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Radius 



Figure 12: The evolution of r 2 p (where p is the energy density of the complex scalar field) with 
massive field (top) and massless (bottom). In the massive case, much of the scalar field collapses 
and a perturbed boson star is formed at the center, settling down by gravitational cooling. In the 
massless case, the scalar field bounces through the origin and then disperses without forming any 
compact object. From [199 . 



4.2 Dynamics of binary boson stars 

The dynamics of binary boson stars is sufficiently complicated that it generally requires numerical 
solutions. The necessary lack of symmetry and the resolution requirement dictated by the harmonic 
time dependence of the scalar field combine so that significant computational resources must be 
expended for such a study. However, boson stars serve as simple proxies for compact objects 
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without the difficulties (shocks and surfaces) associated with perfect fluid stars, and, as such, binary 
BS systems have been studied in the two-body problem of general relativity. When sufficiently 
distant from each other, the precise structure of the star should be irrelevant as suggested by 
Damour's "effacement theorem" [62] . 

First attempts at binary boson star simulations assumed the Newtonian limit, since the SP 
system is simpler than the EKG one. Numerical evolutions of Newtonian binaries showed that in 
head-on collisions with small velocities, the stars merge forming a perturbed star [511 . With larger 
velocities, they demonstrate solitonic behaviour by passing through each other, producing an inter- 
ference pattern during the interaction but roughly retaining their original shapes afterwards [52] , 
In [51] it was also compute preliminary simulation of coalescing binaries, although the lack of res- 
olution in these 3D simulations prevents of strong statements on these results. The head-on case 
was re-visited in [27] with a 2D axisymmetric code. In particular, these evolutions show that final 
state will depend on the total energy of the system, that is, the addition of kinetic, gravitational 
and self-interaction energies. If the total energy is positive, the stars exhibit a solitonic behaviour 
both for identical (see Figure [T3|) as for non-identical stars. When the total energy is negative the 
gravitational force is the main driver of the dynamics of the system. This case produces a true 
collision, forming a single object with large perturbations which slowly decays by gravitational 
cooling, as displayed in Figure [TH 
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Figure 13: Collision of identical boson stars with large kinetic energy in the Newtonian limit. The 
total energy (i.e., the addition of kinetic, gravitational and self-interaction) is positive and the 
collision displays solitonic behaviour. Contrast this with the gravity-dominated collision displayed 
in Figure [Til From [27] . 

The first simulations of boson stars with full general relativity were reported in [13] . where 
the gravitational waves were computed for a head-on collision. The general behaviour is similar 
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Figure 14: Collision of identical boson stars with small kinetic energy in the Newtonian limit. The 
total energy is dominated by the gravitational energy and is therefore negative. The collision leads 
to the formation of a single, gravitationally bound object, oscillating with large perturbations. 
This contrasts with the large kinetic energy case (and therefore positive total energy) displayed in 
Figure [H From [27]. 



to the one displayed for the Newtonian limit; the stars attract each other through gravitational 
interaction and then merge to produce a largely perturbed boson star. However, in this case the 
merger of the binary was promptly followed by collapse to a black hole, an outcome not possible 
when working within Newtonian gravity instead of general relativity. Unfortunately, very little 
detail was given on the dynamics. 

Much more elucidating was work in axisymmetry |143] , in which head-on collisions of identical 
boson stars were studied in the context of critical collapse (discussed in Section 16. ip with general 
relativity. Stars with identical masses of M = 0.47 sa 0.75 M max were chosen, and so it is not 
surprising that for small initial momenta the stars merged together to form an unstable single star 
(i.e., its mass was larger than the maximum allowed mass, A/ max ). The unstable hypermassive 
star subsequently collapsed to a black hole. However, for large initial momentum the stars passed 
through each other, displaying a form of solitonic behaviour since the individual identities were 
recovered after the interaction. The stars showed a particular interference pattern during the 
overlap, much like that displayed in Figures [T] and Q2] 

Another study considered the very high speed, head-on collision of BSs [56]. Beginning with 
two identical boson stars boosted with Lorentz factors ranging as high as 4, the stars generally 
demonstrate solitonic behaviour upon collision, as shown in the insets of Figure [T8l This work is 
further discussed in Section 16.21 

The interaction of non-identical boson stars was studied in [175] using a 3D Cartesian code to 
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simulate head-on collisions of stars initially at rest. It was found that, for a given separation, the 
merger of two stars would produce an unstable star which collapses to a black hole if the initial 
individual mass were M > 0.26 « QAM max . For smaller masses, the resulting star would avoid 
gravitational collapse and its features would strongly depend on the initial configuration. The 
parameterization of the initial data was written as a superposition of the single boson star solution 
c/)q (r), located at different positions ri and r 2 

<t> = (1 Vo(ri)e lwt + (2 Vo(r 2 )e l(ewt+5) . (82) 

Many different initial configurations are possible with this parameterization. The precise solution 
(f>o is unaffected by changing the direction of rotation (within in the complex plane isospace) via 
e = ±1 or by a phase shift S. 

When e = — 1 the Noether charge changes sign and the compact object is then known as an 
anti-boson star. Three particular binary cases were studied in detail: (i) identical boson stars 
(e = 1, S = 0), (ii) the pair in phase opposition (e = 1, 8 = 7r), and (iii) a boson-anti-boson pair 
(e = — 1, S = 0). The trajectories of the centers of the stars are displayed in Figure [151 together 
with a simple estimate of the expected trajectory assuming Newtonian gravity. The figure makes 
clear that the merger depends strongly on the kind of pair considered, that is, on the interaction 
between the scalar fields. 

A simple energy argument is made in |175j to understand the differing behavior. In the weak 
gravity limit when the stars are well separated, one can consider the local energy density between 
the two stars. In addition to the contribution due to each star separately, a remaining term A 
results from the interaction of the two stars and it is precisely this term that will depend on the 
parameters e and S. This term takes the simple form 

A = A cos[(l -e)wi- 6] (83) 

where Ao is a positive definite quantity. One then observes that the identical pair will have an 
increased energy density A = +Ao resulting in a deeper (and more attractive) gravitational well 
between the stars. In contrast, the pair with opposite phases has a decreased energy density 
A = — Ao between them, resulting in a gravitational well less attractive than the area surrounding 
it resulting in repulsion. The boson-anti-boson pair has an interaction that is harmonic in time 
A = Ao cos (2u)t) and therefore sometimes positive and sometimes negative. However, if the time 
scale of interaction is not particularly fast, then the interaction averages to zero. Note that the 
boson-anti-boson pair trajectory is the closest to the simple Newtonian estimate. The qualitative 
behaviour agrees very well with the numerical results. 

The orbital case was later studied in [174] , This case is much more involved both from the 
computational point of view (i.e., there is less symmetry in the problem) and from the theoretical 
point of view, since for the final object to settle into a stationary, rotating boson star it must 
satisfy the additional quantization condition for the angular momentum of Eq. (|71l) . 

One simulation consisted of an identical pair each with individual mass M = 0.5, with small 
orbital angular momentum such that J < N. In this case, the binary merges forming a rotating 
bar that oscillates for some time before ultimately splitting apart. This can be considered as a 
scattered interaction, which could not settle down to a stable boson star unless all the angular 
momentum was radiated. 

In the case of boson-anti-boson pair, the total Noether charge is already trivial, and the final 
object resembles the structure of a rotating dipole. The pair in opposition of phase was not 
considered because of the repulsive effect from the interaction. The cases with very small angular 
momentum J <C N or with J < N collapsed to a black hole soon after the merger. The trajectories 
for this latter case are displayed in Figure [HO indicating that the internal structure of the star is 
irrelevant (as per the effacement theorem [62]) until that the scalar fields overlap. 
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Figure 15: The positiQn of the center of^)QA3S in a head-on QflSy as a function of time for 
(i) [B-B] identical BSs, (ii) [B-poB] opposite phasfe: pair, and (hi) [B-aB] a boson-anti-boson pah. 
A simple argument is made which qualitatively matches these numerical results, as discussed in 
Section 14.21 Also shown is the expected trajectory from a simple Newtonian two-body estimate. 
From [175] . 
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Figure 16: The position of the center of one BS within an orbiting binary as a function of time for 
the two cases: (i) [B-B] identical BSs and (ii) [B-poB] opposite phase pair. Notice that the orbits 
are essentially identical at early times (and large separations), but that they start to deviate from 
each other on closer approach. This is consistent with the internal structure of each member of 
the binary being irrelevant at large separations. From [174J . 
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Other simulations of orbiting, identical binaries have been performed within the conformally flat 
approximation instead of full GR, which neglects gravitational waves (GW) [168] . Three different 
qualitative behaviours were found. For high angular momentum, the stars orbit for comparatively 
long times around each other. For intermediate values, the stars merged and formed a pulsating 
and rotating boson star. For low angular momentum, the merger produces a black hole. No 
evidence was found of the stars splitting apart after the merger. 
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5 Boson Stars in Astronomy 



Scalar fields are often employed by astronomers and cosmologists in their efforts to model the 
Universe. Most models of inflation adopt a scalar field as the inflaton field, the vacuum energy 
of which drives the exponential inflation of the Universe. Dark energy also motivates many scalar 
field models, such as k-essence and phantom energy models. It is therefore not surprising that 
boson stars, as compact configurations of scalar field, are called upon to provide consequences 
similar to those observed. 

5.1 As astrophysical stellar objects 

We have already discussed a number of similarities between boson stars and models of neutron 
stars. Just as one can parameterize models of neutron stars by their central densities, one can 
consider a I-parameter family of boson stars according to the central magnitude of the scalar field. 
Considering the mass as a function of this parameter, one finds the existence of a local maximum 
across which solutions transition from stable to unstable, just as is the case for neutron stars. 
Models of neutron stars can be constructed with different equations of state, whereas boson stars 
are constructed with differing scalar field potentials. 

One difference of consequence concerns the stellar surface. Neutron stars of course have a surface 
at which the fluid density is discontinuous, as discussed for example in |101U103| . In contrast, the 
scalar field that constitutes the boson star is smooth everywhere and lacks a particular surface. 
In its place, one generally defines a radius which encompasses some percentage (e.g. 99%) of the 
stellar mass. Such a difference could have observational consequences when matter accretes onto 
either type of star. 

It is still an open question whether some of the stars already observed and interpreted as neutron 
stars could instead be astrophysical boson stars. In a similar fashion, it is not known whether 
many, if not all, of the stars we observe already have a bosonic component that has settled into 
the gravitational well of the star (see Section [3TB1 for a discussion of fermionic-bosonic stars). The 
bosonic contribution may arise from exotic matter which could appear at high densities inside the 
neutron star or from some sort of dark matter accretion. This possibility has gained popularity 
over the last years and there have been several attempts to constrain the properties of weakly 
interacting dark matter particles (WIMPs) by examining signatures related to their accretion 
and/or annihilation inside stars (for instance, see [139. and works cited in the introduction). 

Recently, it was suggested that, due to the stronger gravitational field of neutron stars compared 
to other stars such as white dwarfs and main sequence stars, WIMPs will accrete more efficiently, 
leading to two different possibilities. If the dark matter is its own antiparticle, it will self-annihilate 
and heat the neutron star. This temperature increase could be observable in old stars, especially 
if they are close to the galactic center |139[ [65] . If WIMPs do not self-annihilate, they will settle 
in the center of the star forming a sort of fermionic-bosonic star. The accretion of dark matter 
would then increase the star's compactness until the star collapses [65] . 

Because of the similarities between boson stars and neutron stars, one finds that boson stars are 
often used in place of the other. This is especially so within numerical work because boson stars 
are easier to evolve than neutron star models. One can, for example compare the gravitational 
wave signature of a boson star merger with that of more conventional compact object binaries 
consisting of BHs and NSs (for a review of BH-NS binaries see |202| ). At early times, the precise 
structure of the stars is irrelevant and the signatures are largely the same. However, for the late 
stages of merger, the relative phase of the boson stars determines the GW signature |175[ 1174] , 

Ref. |172] follows such work by considering the result of a collision between a BH and a boson 
star. In particular, they consider the problem as a perturbation of a black hole via scalar accretion 
and analyze the resulting gravitational wave output. The hope is that observations of gravitational 
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waves that are expected in a few years from ALIGO/ VIRGO will be able to distinguish the type 
of matter accreting onto a BH. 

With the continued advancement in observation, both in the electromagnetic and gravitational 
spectra, perhaps soon we will have evidence for these questions. At the same time, further study 
of boson stars can help identify possible distinguishing observational effects in these bands. One 
example where knowledge is lacking is the interaction between boson stars with a magnetic field. 
Whereas a neutron star can source its own magnetic field and a neutral star can obtain an induced 
charge when moving with respect to a magnetic field, we are aware of no studies of the interaction 
of boson stars with a magnetic field. 

5.2 Compact alternatives to black holes 

As a localized scalar field configuration, a boson star can be constructed as a non-interacting 
compact object, as long as one does not include any explicit coupling to any electromagnetic or 
other fields. In that respect, it resembles a BH, although it lacks a horizon. Can observations 
of purported BHs be fully explained by massive boson stars? See Ref. |184j for a review of such 
observations. 

Neutron stars also lack horizons, but, in contrast to a boson star, have a hard surface. A hard 
surface is important because one would expect accretion onto such a surface to have observable 
consequences. Can a boson star avoid such consequences? Yuan, Narayan and Rees consider the 
the viability of 1OM boson stars as BH candidates in X-ray binaries [224] . They find that accreting 
gas collects not at the surface (which the star lacks), but instead at the center which ultimately 
should lead to Type I X-ray bursts. Because these bursts are not observed, the case against boson 
stars as black hole mimickers is weakened (at least for BH candidates in X-ray binaries). 

Ref. 107 considers a simplified model of accretion and searches for boson star configurations 
which would mimic an accreting black hole. Although they find matches, they find that light 
deflection about a boson star will differ from the BH they mimic because of the lack of a photon 
sphere. Further work, studies the scalar field tails about boson stars and compares them to those 
of BHs |I55j . If indeed a boson star collapses to a BH, then one could hope to observe the QNM 
of the massive scalar field, as described in |116] . 

Some of the strongest evidence for the existence of BHs is found at the center of most galaxies. 
The evidence suggests supermassive objects (of the order of millions of solar mass) occupying a 
small region (of order an astronomical unit) [33 . While definitive evidence for a BH horizon from 
conventional electromagnetic telescopes is perhaps just on the "horizon" [43], there are those who 
argue the viability of supermassive boson stars at galactic centers [215] . There could potentially be 
differences in the (electromagnetic) spectrum between a black hole and a boson star, but there is 
considerable freedom in adjusting the boson star potential to tweak the expected spectrum [105] . 
However, there are stringent constraints on BH alternatives to Sgr A* by the low luminosity in the 
near infrared [44]. In particular, the low luminosity implies a bound on the accretion rate assuming 
a hard surface radiating thermally, and therefore the observational evidence favors a black hole 
because it lacks such a surface. 

However, the observation of gravitational waves from such objects may be able to distinguish 
BHs from BSs [25]. Such a test would occur in the bandwidth for a space-based observatory 
such as the beleaguered LISA mission. Because BSs allow for orbits within what would otherwise 
be a black hole event horizon, geodesies will exhibit extreme pericenter precession resulting in 
potentially distinguishable gravitational radiation [130] . In any case, observations of supermassive 
objects at the centers of galaxies can be used to constrain the scalar field parameters of possible 
mimickers [T5] . 

There are other possible BH mimickers, and a popular recent one is the gravastar [159] . Com- 
mon among all these alternatives, and most significantly, is the lack of an event horizon. Both 
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gravastars and BSs undergo an ergoregion instability for high spin (J/(GM 2 ) > 0.4) [45]. As 
mentioned above for BSs, gravitational waves may similarly be able to distinguish gravastars from 
BHs [I7B] . 

5.3 As origin of dark matter 

Studies of stellar orbits within various galaxies produce rotation curves which indicate galactic 
mass within the radius of the particular orbit. The discovery that these curves remain flat at 
large radius suggests the existence of a large halo of massive, yet dark, matter which holds the 
galaxy together despite its large rotation. However, what precise form of matter could fulfill the 
observational constraints is still very much unclear. Scalar fields are an often used tool in the 
cosmologist's toolkit, but one cannot have a regular, static configuration of scalar field to serve 
as the halo [179] (see [70] as discussed in Section 16.31 for a discussion of rotating boson stars with 
embedded, rotating BH solutions). Instead, galactic scale boson stars are one possible candidate. 

Boson stars can be matched onto the observational constraints for galactic dark matter ha- 
los [1461 1200] , However, multi-state boson stars that superpose various boson star solutions (e.g. 
an unexcited solution with an excited solution) can perhaps find better fits to the constraints [216] . 
Boson stars at the galactic scale may not exhibit general relativistic effects and can be effectively 
considered as Bose-Einstein condensates (BEC) with angular momentum |186j . 

The solitonic nature of boson stars (see Figure [lj lends itself naturally to the wonderful ob- 
servation of dark matter in the Bullet Cluster [147] . Ref. [148] attempts to determine a minimum 
galactic mass from such a match. 

Interestingly Ref. [20] foregoes boson stars and instead looks for quasi-stationary scalar solutions 
about a Schwarzschild black hole that could conceivably survive for cosmological times. Another 
approach is to use scalar fields for both the dark matter halo and the supermassive, central object. 
Ref. [9 a looks for such a match, but finds no suitable solutions. Quite a number of more exotic 
models viably fit within current constraints, including those using Q-balls |72j . 
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6 Boson Stars in Mathematical Relativity 



Although the experimental foundation for the existence of boson stars is completely lacking, on 
the theoretical and mathematical front, boson stars are well studied. Recent work includes a 
mathematical approach in terms of large and small data [89] , followed up by studying singularity 
formation |153] and uniqueness 90, 152 for a certain boson star equation. In Ref. [49], they study 
radial solutions of the semi-relativistic Hartree type equations in terms of global well posedness. 
Ref. [31] demonstrates stationarity of time periodic scalar field solutions. 

Beyond just existence, however, boson stars are often employed mathematically to study dy- 
namics. Here, we concentrate on a few of these topics which have attracted recent interest. 

6.1 Black hole critical behavior 

If one considers some initial distribution of energy and watches it evolve, generally one arrives at 
one of three states. If the energy is sufficiently weak in terms of its gravity, the energy might end 
up dispersing to larger and large distances. However, if the energy is instead quite large, then 
perhaps it will concentrate until a black hole is formed. Or, if the form of the energy supports it, 
some of the energy will condense into a stationary state. 

In his seminal work [53] . Choptuik considers a real, massless scalar field and numerically evolves 
various initial configurations, finding either dispersion or black hole formation. By parameterizing 
these initial configurations, say by the amplitude of an initial pulse p, and by tuning this parameter, 
he was able to study the threshold for black hole formation at which he found fascinating black hole 
critical behavior. In particular, his numerical work suggested that continued tuning could produce 
as small a black hole as one wished. This behavior is analogous to a phase transition in which 
the black hole mass serves as an order parameter. Similar to phase transitions, one can categorize 
two types of transition that distinguish between whether the black hole mass varies continuously 
(Type II) or discontinuously (Type I). For Choptuik's work with a massless field, the transition 
is therefore of Type II because the black hole mass varies from zero continuously to infinitesimal 
values. 

Subsequent work has since established that this critical behavior can be considered as occurring 
in the neighborhood of a separatrix between the basins of attraction of the two end states. For 
p = p* , the system is precisely critical and remains on the (unstable) separatrix. Similarly other 
models find such threshold behavior occurring between a stationary state and black hole formation. 
Critical behavior about stationary solutions necessarily involve black hole formation " turning-on" 
at finite mass, and is therefore categorized as Type I critical behavior. 

The critical surface therefore appears as a co- dimension 1 surface which evolutions increasingly 
approach as one tunes the parameter p. The distance from criticality \p — p* | serves as a measure 
of the extent to which a particular initial configuration has excited the unstable mode that drives 
solutions away from this surface. For Type II critical behavior, the mass of the resulting black hole 
mass scales as a power-law in this distance, whereas for Type I critical behavior, it is the survival 
time of the critical solution which scales as a power-law. See |102] for a recent review. 

We have seen that boson stars represent stationary solutions of Einstein's equations and thus 
one would correctly guess that they may occur within Type I black hole critical behavior. To look 
for such behavior, Hawley and Choptuik jl 1 lj begin their evolutions with boson star solutions 
and then perturb them both dynamically and gravitationally. They therefore included in their 
evolutionary system a distinct, free, massless, real scalar field which couples to the boson star 
purely through its gravity. 

The initial data therefore consisted of a boson star surrounded by a distant, surrounding shell 
of real scalar field parametrized by the amplitude of the shell. For small perturbations, the boson 
star oscillated about an unstable boson star before settling into a low mass, stable solution (see 
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Figure 1X7]) . For large perturbations, the real scalar field serves to compress the initial star and, 
after a period of oscillating about an unstable boson star, the complex field collapses to a black 
hole. By tuning the initial perturbation, they find a longer and longer lived unstable boson star 
which serves as the critical solution (see Figure 1 1 1 1) . The survival time r obeys a power-law in 
terms of the distance from criticality \p — p* \ 



t oc 7 In \p — p* 



(84) 



where 7 is a real constant that depends on the characteristic instability rate of the particular 
unstable boson star approached in the critical regime. 

One can also consider these BSs in axisymmetry in which non-spherically symmetric modes 
could potentially become important. A first step in this direction studied spherically symmetric BSs 
within conformally flat gravity (which does not allow for gravitational waves) in axisymmetry 1 88 . 
Later, better resolution using adaptive mesh refinement within full general relativity was achieved 
by |143[ 1144] which upheld the results found within spherical symmetry. This work thus suggests 
that there are either no additional, unstable, axisymmetric modes or that such unstable modes are 
extremely slowly growing. 
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Figure 17: Evolution of a boson star (solid line) perturbed by a shell of scalar field (dashed line). 
Shown is the mass density dM/dr for each contribution. By t w 100 the real scalar field pulse has 
departed the central region and perturbed the boson star into an unstable, compact configuration. 
Contrast the t = frame with that of t = 97.5 and note the increase in compaction. This unstable 
BS survives until t rj 500 only because the initial perturbation has been tuned to one part in 10 15 
and indicates Type I critical behavior. From |144j . 

A very different type of critical behavior was also investigated by Lai [143] , By boosting 
identical boson stars toward each other and adjusting their initial momenta, he was able to tune 
to the threshold for black hole formation. At the threshold, he found that the time till black hole 
formation scaled consistent with Type I critical behavior and conjectured that the critical solution 
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was itself an unstable boson star. This is one of the few fully nonlinear critical searches in less 
symmetry than spherical symmetry, and the first of Type I behavior in less symmetry. A related 
study colliding neutron stars instead of boson stars similarly hnds Type I critical behavior [126 
and subsequently confirmed by [129 . 

The gauged stars discussed in Section 13.91 also serve as critical solutions in spherical symme- 
try pi 1551 1167]. 

6.2 Hoop conjecture 

An interesting use of boson stars was made by Choptuik and Pretorius [56] . They sought to answer 
classically whether the ultra-relativistic collision of two particles results in black hole formation. 
Such a question clearly has relevance to hopes of producing black holes at the LHC (see for 
example |145| 1178] ). Guidance on this question is provided by Thome's Hoop Conjecture [214 
which suggests that, if one squeezes energy into some spherical space of dimension less than the 
Schwarzschild radius for that energy, then a black hole is formed. 
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Figure 18: Evolutions of the head-on collisions of identical boson stars boosted toward each other 
with initial Lorentz factors 7 as indicated. Time flows downward within each column and the top 
edge displays the axis of symmetry. The color-scale indicates the value of \<f>\. In the middle frames 
one sees the interference pattern characteristic of high kinetic energy BS collisions (as mentioned 
in Figure [1]). In the last column on the right, the collision produces a BH with apparent horizon 
indicated by the black oval in the third frame. From |56| . 

They therefore numerically collide boson stars head-on at relativistic energies to study black 
hole formation from just such dynamical "squeezing" . Here, the nature of boson stars is largely 
irrelevant as they serve as simple bundles of energy that can be accelerated (see Figure IT8|) . How- 
ever, unlike using boosted black hole solutions, the choice of boson stars avoids any type of bias or 
predisposition to formation of a black hole. In addition, a number of previous studies of boson star 
head-on collisions showed interesting interference effects at energies below the threshold for black 
hole formation [501 1521 [1431 1168J . Indeed, it has been proposed that such an interference pattern 
could be evidence for the bosonic nature of dark matter because of evidence that an ideal fluid 
fails to produce such a pattern [97] . 

Choptuik and Pretorius |56j find that indeed black hole formation occurs at energies below 
that estimated by the Hoop Conjecture. This result is only a classical result consistent with the 
conjecture, but if it had not held, then there would have been no reason to expect a quantum 
theory to be consistent with it. 
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6.3 Other dimensions 



Much work has been invested recently in considering physics in other dimensions. Motivation comes 
from various ideas including string theory (more dimensions) such as the AdS/CFT correspondence 
and holography (one fewer dimensions) [1581 11601 1182] . Another source of motivation comes from 
the fact that higher dimensional black holes can have very different properties than those in three 
spatial dimensions [75]. Perhaps BSs will similarly display novel properties in other dimensions. 

In lower dimensional AdS (2 + 1) spacetimes, early work in 1998 studied exact solutions of boson 
stars 192, 68] 1193] . Higher dimensional scenarios were apparently first considered qualitatively a 
few years later in the context of brane world models [206]. This discussion was followed with a 
detailed analysis of the 3, 4, and 5 dimensional AdS solutions |12) . 

More recently, Ref. [81] considers oscillatons in higher dimensions and measure the scalar mass 
loss rate for dimensions 3, 4, and 5. They extend this work considering inflationary space-times 80 . 

The axisymmetric rotating BSs discussed in Section ^. 5l satisfv a coupled set of nonlinear, elliptic 
PDEs in two dimensions. One might therefore suspect that adding another dimensions will only 
make things more difficult. As it turns out, however, moving to four spatial dimensions provides 
for another angular momentum, independent of the one along the z-direction (for example). Each 
of these angular momenta are associated with their own orthogonal plane of rotation. And so if 
one chooses solutions with equal magnitudes for each of these momenta, the solutions depend on 
only a single radial coordinate. This choice results in the remarkable simplification that one need 
only solve ODEs to find rotating solutions [140] , 

In [110J . they extend this idea by assuming an ansatz for two complex scalar fields with equal 
magnitudes of angular momentum in the two independent directions. Letting the complex doublet 
be denoted by the ansatz takes the form 



in terms of the two angular coordinates tpi and tp<i. One observes that the BS (i) retains a profile 
<p(r), (ii) possesses harmonic time dependence, and (iii) maintains single- valuedness in the two 
angles (the ansatz assumes a rotational quantum number of one). They find solutions that are 
both globally regular and asymptotically flat but these solutions appear only stable with weak 
gravitational coupling [110] . 

The work of [70] makes ingenious use of this 5D ansatz to construct rotating black holes with 
only a single killing vector. They set the potential of [110] to zero so that the scalar fields are 
massless and they add a (negative) cosmological constant to work in anti-de Sitter (AdS). They 
find solutions for rotating black holes in 5D AdS that correspond to a bar mode for rotating neutron 
stars in 3D (see also [203 for a numerical evolution of a black-hole in higher dimensions which 
demonstrates such bar formation). One might expect such a non-symmetric black hole to settle 
into a more symmetric state via the emission of gravitational waves. However, AdS provides for 
an essentially reflecting boundary in which the black hole can be in equilibrium. The distortion of 
the higher dimensional black hole also has a correspondence with the discrete values of the angular 
momentum of the corresponding boson star. For higher values of the rotational quantum number, 
the black hole develops multiple "lobes" about its center. 

These solutions are extended to arbitrary odd-dimensional AdS spacetimes in 207 . Finding 
the solutions perturbatively, they explicitly show that these solutions approach (i) the boson star 
and (ii) the Myers-Perry black hole solutions in AdS 171 in different limits. See [76] for a review 
of black holes in higher dimensions. 

The same authors of 70^ also report on the existence of geons in 3+1 AdS "which can be 
viewed as gravitational analogs of boson stars" [71] (recall that boson stars themselves arose from 
Wheeler's desire to construct local electrovacuum solutions). These bundles of gravitational energy 
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are stable to first order due to the confining boundary condition adopted with AdS. However, these 
geons and the black hole solutions above [70] are unstable at higher order because of the turbulent 
instability reported in [32] . 

Ref. [22] also studies black hole solutions in 5D AdS. They find solutions for black holes with 
scalar hair that resemble a boson star with a BH in its center. See Ref. [94] for a review of charged 
scalar solitons in AdS. 

7 Final remarks 

Boson stars have a long history as candidates for all manner of phenomena, from fundamental 
particle, to galactic dark matter. A huge variety of solutions have been found and their dynamics 
studied. Mathematically, BS are fascinating soliton-like solutions. Astrophysically, they represent 
possible explanations of black hole candidates and dark matter, with observations constraining BS 
properties. 

Further constraints appear to be right around the corner in a few directions. The LHC is 
rapidly narrowing the possible mass range for a Higgs particle (the quantized version of a funda- 
mental scalar field), and some hold hope that black holes might be produced that would indicate 
the existence of higher dimensions. COGENT, Pamela, and others are finding intriguing (and frus- 
tratingly contradictory) clues about the true nature of dark matter. And Planck and JWST, if it 
can overcome its funding and budgetary problems, promise to further refine our view of the cosmos 
and the role of scalar fields within. Advanced LIGO will be complete in a few years and the future 
of space-based GW telescopes such as LISA is currently being determined. GW observations will 
be ground-breaking on so many fronts, but could potentially help distinguish BSs from neutron 
stars or black holes. 

Perhaps future work on boson stars will be experimental, if fundamental scalar fields are ob- 
served, or if evidence arises indicating the boson stars uniquely fit galactic dark matter. But 
regardless of any experimental results found by these remarkable experiments, there will always 
be regimes unexplored by experiments where boson stars will find a natural home. 
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